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We investigate the non-linear coupling between radial and non-radial oscillations of static spheri- 
cally symmetric neutron stars as a possible mechanism for the generation of gravitational waves that 
may lead to observable signatures. In this paper we concentrate on the axial sector of the non-radial 
perturbations. By using a multi-parameter perturbative framework we introduce a complete descrip- 
tion of the non-linear coupling between radial and axial non-radial oscillations; we study the gauge 
invariant character of the associated perturbative variables and develop a computational scheme to 
evolve the non-linear coupling perturbations in the time domain. We present results of simulations 
corresponding to different physical situations and discuss the dynamical behaviour of this non-linear 
coupling. Of particular interest is the occurrence of signal amplifications in the form of resonance 
phenomena when a frequency associated with the radial pulsations is close to a frequency associated 
with one of the axial w-modes of the star. Finally, we mention possible extensions of this work and 
improvements towards more astrophysically motivated scenarios. 
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I. INTRODUCTION 

Gravitational wave astronomy is currently developing 
as a new way of studying the cosmos, with a num- 
ber of different experiments starting operating or be- 
ing developed: earth-based interferometers (LIGO 0, 
VIRGO 0, GEO600 and TAMA ||), resonant bars 
(EXPLORER, AURIGA, NAUTILUS, ALLEGRO and 
NIOBE; see @) and spheres such as MiniGRAIL as 
well as the laser interferometer space antenna LISA 0. 
The scientific success of these detectors in producing new 
physical and astrophysical knowledge depends on the 
amount of a priori available theoretical understanding 
about the different sources of gravitational waves. It is 
then of outmost importance to have good theoretical de- 
scriptions of the systems that are likely to be observed 
by the ongoing gravitational- wave experiments. 

Due to the weak character of the gravitational inter- 
actions there is only a few number of systems that can 
generate gravitational waves detectable for the experi- 
ments just mentioned. In general, the dynamics of these 
systems involves strong gravitational fields that need to 
be described, using different levels of approximations, by 
general relativity. This is the case in the modeling of com- 
pact objects such as neutron stars and supernova? core 
collapse, where non-linearity plays a role in the dynam- 
ics. For these systems, relativistic perturbation theory is 



an excellent approach. The linear theory has been used 
for a long time to study their oscillations and instabili- 
ties However, we know relatively little about non- 
linear dynamical effects (mainly due to numerical stud- 
ies [lSlillllEQISElIa). Despite the fact that 
strong non-linear effects require a fully non-perturbative 
approach, it is reasonable to expect that there are in- 
teresting physical phenomena that only involve a mild 
non-linearity for which a second-order treatment should 
be perfectly suited. In this sense, a non-linear pertur- 
bative approach is therefore timely and may lead to a 
better understanding of the mechanisms of generation of 
gravitational waves. 

An interesting scenario to study is a neutron star which 
is oscillating radially and non-radially. At first order 
radial oscillations of a spherical star don't emit per se 
any gravitational waves, but they can drive and possibly 
amplify non-radial oscillations and then produce gravita- 
tional radiation to a significant level. In addition one may 
expect the appearance of non-linear harmonics, which 
may also come out at lower frequencies than the linear 
modes ^3 ^| , where the Earth-based laser interferom- 
eters have a higher sensitivity. An additional motivation 
is that there are a number of studies aiming at investi- 
gating if non-radial oscillations of stars can be excited 
by external sources (see e.g. 0,H3). However, our idea 
is to investigate whether non-radial oscillations can be 
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driven or even amplified through coupling by an inter- 
nal radial oscillation, regardless of the presence of an 
external source. These non-linear processes can occur, 
for instance, in a proto-neutron star that is still pulsat- 
ing. A mainly radial pulsation could for example drive 
the non-radial oscillations, cither naturally present, or 
excited through fall-back accretion. 

In recent papers |2ll |22( a non-linear multi-parameter 
perturbative formalism has been developed, which can 
be applied to the study of the non-linear dynamics of the 
oscillations of neutron stars, as shown in j2^| where it 
has been further developed to study a particular second- 
order effect: the coupling of radial and non-radial first 
order perturbations of a compact spherical star. One 
of the main ideas behind this work is that it is very 
convenient to treat the two sets of perturbations, ra- 
dial and non-radial, as separately parametrized and then 
to study, at the next perturbative order, the way in 
which they couple, which corresponds to a particular sec- 
tor of the second-order perturbations, the one formed 
by the product of radial and non-radial perturbations. 
In |23|| . the equations describing the coupling were ob- 
tained for the case of polar non-radial perturbations and 
gauge-invariant variables were also found (having fixed 
the gauge for the radial ones for simplicity). 

In this paper, we extend this study to the case of axial 
perturbations. At first order, axial perturbations are de- 
coupled from fluid perturbations, but at the next pertur- 
bative order, the terms describing the coupling are driven 
by the radial pulsations. Here we discuss the results of 
simulations in the time domain that show how the cou- 
pling mechanism works and under which conditions we 
may have a situation in which this coupling produces 
amplifications due to resonant phenomena. 

The plan of this paper is the following: In Section ITT1 
we describe the perturbative framework we use in this 
work. In Section IIIII we introduce all the necessary in- 
gredients required by the perturbative scheme in order 
to study the coupling of radial and axial non-radial os- 
cillations of neutron stars. In Section Hvl we describe the 
structure of the computational framework we have imple- 
mented to solve the perturbative equations in the time 
domain. In particular, we discuss in detail the construc- 
tion of initial data and the numerical treatment of the 
boundary conditions. In Section we present the re- 
sults of our simulations for two different physical scenar- 
ios: (i) the scattering of a gravitational wave packet by 
a differentially rotating star, and (ii) the non- linear cou- 
pling in a differentially rotating and radially oscillating 
compact star. In particular we discuss the appearance 
of amplifications due to resonances of the system. We 
conclude in Section IVII where we summarize the main 
results of this paper and discuss potential extensions of 
the approach we have used. In the Appendices we give 
some proofs of the gauge invariance character of some 
of the perturbative quantities (Appendix^), expressions 
for the source terms in the perturbative equations de- 
scribing the coupling (Appendix [BJ) , and the equations 



for radial perturbations in terms of the tortoise fluid co- 
ordinate (Appendix Q . 

Throughout this work we use the following conven- 
tions: We use Greek letters to denote spacetime indices; 
capital Latin letters for indices in the time-radial part 
of the metric; lower-case Latin indices for the spheri- 
cal sector of the metric. We use physical units in which 
8ttG = c = 1 . 



II. PERTURBATIVE FRAMEWORK 

In order to investigate the effects of the coupling of 
linear axial perturbations and first order radial oscilla- 
tions we use a general two-parameter framework recently 
developed [2lH22| . In essence, this consists of separately 
parametrizing the radial and non radial oscillations. We 
have presented the formalism for the coupling between 
radial and polar non-radial perturbations in |23j . In this 
paper we develop and apply the formalism to study the 
coupling arising from axial perturbations, following the 
same general lines. In the following we recall the general 
two-parameter framework. However, to fix the ideas, we 
will refer explicitly to its application to the case of stellar 
oscillations. 

The physical spacetime of the oscillating star we con- 
sider is represented as a non-linear perturbation of a 
static spherical background stellar model, with metric 
g(o,o) rpj^ ex p ans i on f th e physical metric g around 
this background has the following form: 

„ _ 0.(0,0) , (1,0) (0,1) - (U), n A2 J2s m 

where A and e are the expansion parameters associated 
with radial and non-radial perturbations respectively. 
We use superscripts (I, J) to denote perturbations of or- 
der A / e J . Therefore, the terms g^ 1 ' ) and g^- 1 ' in (JTJ 
represent, respectively, first-order radial and non-radial 
perturbations, while g^ 1,1 ' is the non-linear contribution 
due to their coupling. This is the quantity we are inter- 
ested in this paper, while we will neglect the self-coupling 
terms of order A 2 and e 2 . 

Any physical quantity can be expanded as the metric 
in ft)). In particular, the energy-momentum tensor has 
the expansion: 

Tap = T^+\T^ + eT^ + \eT^ + 0(\*,e 2 ) . (2) 

Let us now analyze the structure of the field equations 

E[g,^ A ]=G[g]-T[g,^ A } = (3) 

arising from the perturbative expansion (see e.g. |24j|L 
Here G denotes the Einstein tensor, T the energy- 
momentum tensor and if) A (A=l, . . . ) the matter vari- 
ables. Introducing (|TJ and @ into Eq. we obtain the 
perturbative equations at each order: 
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E[g,ip A ] = E 
+ XeE^ 



g ( °' 0) ,^ 



; (M) jV ,(M) 



(0,0) 
(1,0) 



A£ 



(1,0) 



g 



e ^(o,i) 



,(0,1) 



(0,1) 



A) 



,(0,1) 



+ O(A 2 ,e 2 )=0. (4) 



This equation has to be satisfied for arbitrary values of 
the expansion parameters, therefore each perturbative 
term must vanish independently. This naturally gives 
rise to an iterative scheme, where at each non-linear or- 
der one has to solve inhomogeneous equations with source 
terms containing the solutions of the lower order equa- 
tions. 

The lowest order in ®, E (Q - 0) = 0, gives the equa- 
tions for the background, in our case the Tolman- 
Oppcnhcimcr-Volkoff (TOV) equations. The other terms 
fli 1 !-!) _ o re p resen t the perturbative equations of or- 
der X 1 e J . The structure of the operators E^' J \ with 
(/, J) (0, 0), is such that they act linearly on every term 
specified between the square brackets, while they are non- 
linear functions of the background quantities g( 0,0 ) and 

we obtain the equations de- 



At first order in A, 



scribing the radial perturbations of the TOV background, 
£;(i,o) _ wri ereas axial non-radial perturbations come 
from the first-order terms in e: 



^(0,1) 



,(o,i) 



0. 



(5) 



Finally, the equations of order Ae describe the radial non- 
radial coupling: 



E 



(i,i) 



,(i,i) 



g 



(1,0) 



g 



,(0,1) 



0, 



(6) 



where the vertical bar separates the unknowns of order 
Ae from the quantities computed in the previous iteration 
of the perturbative scheme, which form the source terms 
of the equations. 

In order to solve these equations, it is important 
to take into account that the operator E^- 1 ' 1 ' acts on 
(g(M) ^U' 1 ') i n Eq. (JSJ in the same way as f^ ' 1 ) acts 
on (gC 3 * 1 ' , i/j^' 1 -*) in Eq. (JSJ). The reason is that both op- 
erators come from the linearization of the Einstein tensor 
operator acting on axial non-radial perturbations. Then, 
using the linearity of E^ 1 ' 1 ^, we can define 



£(bi)[.|o] =J E(°.i)[.] 



as the axial non-radial perturbation operator, 
equation JBJ can be written in the form: 



(7) 
Hence, 



'NR 



,(1,1) )V ,(V) 



(1,0) 



,(0,1) 



,g 



,(0,1) 



(8) 



This particular structure of the equations is very im- 
portant in solving them by using time-domain numerical 
methods. Indeed, we can easily extend a numerical code 
that solves the equations for first-order axial non-radial 
perturbations simply by adding the source terms on the 
right-hand side of (JSJ) . Later on we will discuss additional 
advantages of this structure. 



III. PERTURBATIVE SCHEME FOR STELLAR 
OSCILLATIONS 

In this section we introduce the basic ingredients 
needed to study the coupling of radial and axial per- 
turbations of relativistic stars: the background descrip- 
tion fSec. llTTA^I . the first-order radial and axial pertur- 
bations fSecs. UnBl and lill C|l . and finally the equations 
describing the coupling of the first-order perturbations 
(Sec. lrTlD|) . To simplify our expressions we denote the 
background quantities with an overbar, e.g. g a0 = 



and T 



a/3 



T, 



(0,0) 



a/3 



A. The static spherically-symmetric background 

The background is an equilibrium configuration de- 
scribed by a static spherically-symmetric metric (see, 
e.g., a TOV model. The line-element is: 



Ea(3 



dx a dx 13 



. 2* 



z 2A dr 2 



r 2 (d9 2 



where $ and A are functions of r only. This is a solu- 
tion of Einstein's equations with a perfect-fluid energy- 
momentum tensor 

T a p = (p+p)u a U +pg a0 , (10) 

where the four-velocity of the fluid is u a = (— e , 0, 0, 0) , 
and in general p and p denote the energy density and the 
pressure. The mass function m(r) is introduced through 
the metric function A: e ~ 2A M = l-2m(r)/r. The TOV 
equations of stellar equilibrium read: 



4> 



47T7 • p 



r(r — 2m) 



P.r 

P+P 



Airr p . 



(11) 
(12) 



Specifying the Equation of State (EoS) of the stellar 
equilibrium configuration yields a 1-parameter family 
of solutions of the equations (|11I12|I . which is typically 
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parametrized by the value of the central density p c = 
p(r = 0) . Here we consider a configuration described by 
an adiabatic, barotropic EoS generalization of the New- 
tonian polytropic EoS: 

P = Kf , (13) 

where K is a constant and T the polytropic index. In 
particular, we take K = 100 km 2 , r = 2 , p c = 3 x 10 15 
g/cm 3 . These values yield a stellar configuration with 
mass M = 1.26M Q and radius R = 8.86 km. We denote 
by c s = (dp/dp) 1 / 2 the speed of sound. 

B. Radial Perturbations of the TOV Model 

Radial perturbations (first investigated by S. Chan- 
drasekhar j2|| ; see also ) , by definition, preserve 
the spherical symmetry of the TOV background. They 
can be completely described in terms of the Lagrangian 
displacement of a fluid element around its equilibrium 
position. Alternatively, in the radial gauge, they can be 
written in terms of two metric functions: 

g^ 0) dx a dxP = e 8 *^ 1 ' ) - 2 V W)dt 2 

+ re 2A S^dr 2 . (14) 



The radial perturbations of the fluid variables in this 
gauge are: 

u^dx a = e*(§5( 1 '°)-»?( 1 ' ))A + e A 7( 1 ' )dr,(15) 
5pW = pu>W , (16) 
SpW = c 2 s puj^ a K (17) 



To simplify the equations we use the following pertur- 
bative quantities: metric variables: the metric perturba- 
tions rf 1 '®} and S^ 1,0 ^ ; matter variables: the fluid velocity 
perturbation -y(i>°) and the enthalpy perturbation, which 
is related to the density perturbations through the rela- 
tion: 



#(i,o) = _3p_ w (i,o) . (i 8 ) 
P + P 



The dynamics of radial pulsations is then governed by 
the following system of PDEs: 



< 0) 



7^ 



-cie*" A \ 7 ( ^ 0) + e 2A 



4nrp 



-2A 



— 47r(p + p) 



„(b0) 



-87r(p + p)e* +A 7 (1 ^ 0) , 

I 



(19) 

(20) 
(21) 



which was derived in using the formalism of 
El EM EI3 (see Section UlI CI below^l. A similar system 
was obtained in [34| using the well-known ADM 3+1 for- 
malism [3^, |3(j. Having solved the above system, the 
metric variable jjl 1 ' ) can then be obtained from the fol- 
lowing relation 



rj^' 0) = 4Trr(p+p) 



,5(1.0) 



1 



1 



#(i,o) 



2 A 



(22) 

Finally, our variables have to satisfy the Hamiltonian 
constraint, which takes the following form 



0(1,0) 

w r 



„2A 



I 8irrp - 



2 m 



(23) 

Boundary conditions must be fixed at the center and 
at the stellar surface r = R. The vanishing at the sur- 
face of the first-order Lagrangian pressure perturbation, 
Ap( 1,0 )(i?) = 0, leads to the following condition on ^f 1 ' ) 



(p+p)c 2 s (r 2 e-W m )Jr =R = 0. 



(24) 



The behaviour of S^ 1 ' ) and on the surface can 

be derived from the general evolution equations l|19|) and 
(|21|l . In the radial gauge there is still a residual gauge 
degree of freedom, which can be used for setting to zero 
all the pure gauge perturbations of the vacuum external 
spacetime [33| . For stellar models where the density van- 
ishes on the stellar surface, this gauge choice induces also 
5(1,0) a nd ^(1,0) i Q van ish on that surface |33j|. At r = 0, 
regularity conditions imply 



5 d,o) 
^(1,0) 
#(1,0) 
7 (i,o) 



rS^°\t) + 0(r 3 ) , 
^>°)(i) + 0(r 2 ), 
H^(t) + 0(r 2 ), 
r 7 ( 1 '°)(t)+0(r 3 ). 



(25) 
(26) 
(27) 
(28) 



5 



C. Axial Perturbations of the TOV Model 

Axial oscillations of a static star were first investigated 
by Thorne and Campolattaro [^3, EE > an d have been 
the subject of many subsequent studies (see, e.g. |40ll4l| ). 

Non-radial perturbations are nicely described in the 
formalism introduced by Gerlach and Sengupta [2^, 
l3f1 | and further developed by Gundlach and Martfii- 

Garcfa jumm (which we denote as GSGM formal- 
ism). This was originally introduced to study first or- 
der perturbations of a general time-dependent spherically 
symmetric stellar background. In particular we can ap- 
ply it to our case, where the background is the static 
TOV metric. 



where Y lm are the scalar spherical harmonics, in the fol- 
lowing form 



(o,i) 

y a /3 



,(0,1) _ 



l q/3 





h l r s i a r 



St 1 ^ S l a m 



h lm S lrn j 

5t l X> S l a m 

St lm S lm 



(33) 
(34) 



where h l A h lm , St 1 ™ , and 8t lm are functions of t and r 
only. To simplify the notation, in the following we drop 
the harmonic indices (I, m). As shown in [^Ijfflj ]. a com- 
plete set of axial gauge-invariant quantities is obtained 
by taking the following combinations of h A , h , 8t A , and 
St: 



1. Summary of the GSGM formalism for axial 
perturbations 

The key point in the GSGM formalism is the fact 
that the background manifold is a warped product M 2 x 
S 2 , where S 2 denotes the 2-sphere and M 2 is a two- 
dimensional Lorentzian manifold. The metric can be 
written as the semidirect product of a general Lorentzian 
metric on M 2 , g AB , and the unit curvature metric 7 ab on 
S 2 : 



2hv 



A ' 



St, 



p h 



A - 



L 



St — ph , 



(35) 
(36) 
(37) 



where k A and L A are defined for I > 1 , L is defined 
for I > 2, and v A — r~ 1 r^ A . The only axial fluid per- 
turbation function arises from the expansion of the fluid 
velocity 



4°- 1 ) = (0, J 9S o ) ) 



(38) 



= ( 9ab 

icfi \ o r 2 lah 



(29) 



Hereafter x A denote the coordinates on M 2 , x a the co- 
ordinates on S 2 and r = r(x A ) is a function on M 2 that 
coincides with the invariantly defined radial (area) co- 
ordinate of spherically-symmetric spacetime. A vertical 
bar is used to denote the covariant derivative on M 2 and 
a semicolon to denote that on S 2 . The metrics are co- 
variantly conserved, i.e. 9 AB \ C = l a b-c = • O ne can 
introduce the completely antisymmetric covariant unit 
tensors on M 2 and on S 2 , Z AB and e ab respectively. 

The energy-momentum tensor has the same block di- 
agonal structure as the metric: 



T aj3 = diag (pu A u B +pn A n B , r 2 pj ab ) , 



(30) 



where u A are the non-zero components of the fluid veloc- 
ity, u a = (u A , 0) , and n A is a unit spacelike vector given 
by 



where f3 is a gauge invariant variable. As it has been 
shown in j^, for a perfect fluid the energy-momentum 
gauge-invariant quantities take the form 



L A = (p + P)Pu 



A i 



L = 0. 



(39) 



Then, the equations for the perturbations can be decou- 
pled in terms of the following quantities j2£j 



* = r 3 e AB (r- 2 k A ) ]B , $=Q> + p)P 



and their expressions are 



u A f3\ A + (fi + 2U)$ = , 
where the background quantities p. and U are 



(40) 



A)\B ■ 



(41) 
(42) 



-£ab u 



B 



0. 



(31) 



Here we will focus on the axial non-radial perturba- 
tions (see [23| for the polar ones). The metric and energy- 
momentum axial perturbations can be expanded in terms 
of the vector and tensor axial (or odd-parity) spherical 
harmonics 



S. 



lm 



b\rlm 
1 1 -h 



elm 



5lm 
n-h 



cilm 



b:a 1 



(32) 



\A- 



u 



u V j 



— 1 — A 

r u r, 



(43) 



We can solve these equations by prescribing initial data 
for $ and $ at a time t Q , that is, giving (\& 0) * , j3 ) . 
Equation l|42|) is an equation for (3 only that can be solved 
independently from l|41l) . Its solution is constant along 
the integral curves of u A (see |23). With the solution of 
(|42|l we can solve the equation (|41|) for , where (3 enters 
only in the source term. Once we have solved equations 
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(I41I42|I for \& and /3, we can get the gauge-invariant met- 
ric perturbations k A by means of the following relation: 

(I + 2)k A = l^r 2 pu A - e AB (r^ B . (44) 

We emphasize that the formalism as described so far ap- 
plies to a general time dependent, spherically symmetric 
background spacetime. 



2. GSGM formalism on the TOV background 

In the case of our TOV static background, we have: 

u A = (e-*,0), n A = (0,e- A ), p = U = 0. (45) 

Then, the equations for the gauge-invariant axial pertur- 
bations \& and (3 are 



, 2* 



4?r(p - p) 



6M 



1(1 + 1) 



c ZA 3, r - Uvrp- — | i 



M\ 
r~2J 



0, 



(46) 
(47) 



where, to simplify the equations, we have introduced the 
tortoise coordinate 

dr, = e*- A dr . (48) 

Equation l|47|) tells us that /3 is not a dynamical variable 
and hence, at first order, is t ypically of no interest and 
can be consistently set to zero |32i|4(j,|41|. In this sense, it 
is said that first-order axial perturbations do not depend 
on the fluid motion. 

As we discuss later in Section TlV Bl (3 lm can however 
be taken to represent the harmonic component Im of a 
first order differential rotation profile, and as such plays 
a very interesting role in the coupling, at order Ae, with 
radial oscillations. 

Then, if we consider a non-vanishing (3 , the general 
solution of l|46|) is the combination of (i) the solutions of 
the homogeneous associate equation that describes the 
only dynamical degrees of freedom of axial perturbations, 
i.e. the gravitational waves, and (ii) a particular (static) 
solution that is related to the dragging of inertial frames 
induced by a non- vanishing /3, as discussed in |37|. 

Outside the star, equation 1461) reduces to the well- 
known Regge- Wheeler equation 

-¥ itt + * ir . r .-V,(r)¥ = 0. (49) 

where 

™-(i-^)(!S£i-!?), (so, 

and M is the gravitational mass of the TOV star. In (|49H . 
r„ is the usual Regge- Wheeler tortoise coordinate r„ = 
r + 2M\n[r/(2M) — 1] . To complete the description of 
the first-order axial perturbations we need to discuss the 
boundary conditions at the origin, at the stellar surface, 
and at infinity. The requirement of regularity of the per- 
turbations at the origin yields 

f3 „ r l+1 , * - r l+1 . (51) 



The matching conditions on the stellar surface imply the 
continuity of metric variable ^ , of its time derivative, 
and of the following quantity |33(: 

e- A (r^ 3 *) r - 167rr~ 2 /3. (52) 

In the case of a barotropic equation of state, the pressure 
and mass-energy density vanish on the surface, and the 
condition (|52|) reduces to the continuity of \P r . At infin- 
ity we can simply use Sommerfeld outgoing wave condi- 
tions. 

The master function ^ is related to the emitted power 
in gravitational waves at infinity (see, e.g. H^): 

where we have explicitly restored the harmonic indices 
(l,m) and the overdot denotes time differentiation with 
respect to the Schwarzschild coordinate time. 

D. Non-linear Coupling between radial and 
non-radial perturbations 

To derive the equations that describe the non-linear 
coupling between first-order radial and axial non-radial 
oscillations we are going to follow the procedure used 
in the polar case |23|. In brief, the idea is to obtain 
the equations for the coupling from the equations for the 
axial perturbations f4*T)l and , by using the following 
formal perturbative expansion: 

ga/3=gi°/3+ £ gS> ( 54 ) 

where gfp is a time dependent background composed of 
our static TOV background with the radial oscillations, 
that is g^°j = g a0 + Ag^ 0) , and g£j represents axial 
non-radial perturbations containing the first-order radial 



7 



perturbations and the non-linear coupling terms, that is 

§a/3 = So/3 + ^&a/3 • A s we have already mentioned in 
Sec. [n] this is due to the fact that the first-order axial 
and the coupling perturbations have essentially the same 
structure and are governed by the same type of equations. 

The metric and energy-momentum axial perturbation 
of order (1, 1) can be expanded in odd-parity spherical 
harmonics as in the first-order case [see equations (I33L 
1)340] . To distinguish them we use superscripts (I, J) de- 
noting the perturbative order X 1 e J . In this way, the 
(0,1) metric and energy- momentum perturbations are 
denoted by hf 1] , , tf 1] , and ; and the 

(1, 1) ones by h^ 1] , , t^ A) , and . Likewise, 

the velocity perturbations of order (1,1) are written as 

Then, the gauge-invariant quantities that we construct 



by using the ansatz of l|54(l are: 

= h^-hff+2h^v A , (55) 

= St^-ph^-pc^h^, (56) 
L* 1 ' 1 ) = StW - ph^V - pguMhP'V . (57) 

Comparing expressions H55(l - (|57|l with expressions 135(1 - 
(|37[1 we notice that they differ because of the presence of 
the products of first-order perturbations in the (1,1) per- 
turbations in H56I57JI (see also 23]). The gauge invariance 
of the (1, 1) quantities l|55|l - 157|) is shown in Appendix 1X1 

We can also construct the coupling perturbations fit 1 ' 1 ) 
and vlK 1 ' 1 ) in the same way as /3 (oa) = (3 and ^(o. 1 ) = v]/ 
in 140JI (see Appendix lAl for more details). We find that 
they satisfy the following equations: 



ft* 



(i,i) 



(i,i) 



+ e 



2<i> 



' . . 6M 1(1 + 1) 
4n{p -p) + — K —^- 



9<- l 'V + 167rre 2 * H 



2A / 3(l,l) 



Airrp - 



M 



0(1,1} 



(58) 
(59) 



where S^, and are source terms made out of products 
of radial and axial perturbations. Their expressions can 
be found in Appendix [5] In the stellar exterior we do 
not have fluid perturbations and the radial perturbations 
vanish due to Birkhoff's theorem. Hence, the equations 
for the exterior do not have source terms and then 158|) 
reduces to the Regge- Wheeler equation for $( 1 > 1 ) [same 
as for ^'(O' 1 ), see Eqs. (|49I50[) ]. Therefore, the emitted 
power in gravitational waves at infinity due to the (1,1) 
perturbations is given by equation (J53J, substituting $> lm 
by 'j';™ ■ Regarding the boundary conditions for /SK 1,1 ) 
and v]/(i,i) i at the origin they are the same as the ones 
for Z^ - 1 ' and W^' 1 ) [see Eq. l|5T)l]: at infinity we use 
outgoing wave Sommerfeld conditions (see below) ; at the 
stellar surface, both vIa 1 ' 1 ) and its time derivatives are 
continuous. The matching conditions |33j also imply the 
continuity of the following quantity 



+ £^(o.i) _ £*(o,i) 5(1,0) 



2 2 

+ e-V 1 ' 0) * ( *°' 1) + 16™73( 1 < 1 ) 



(60) 



This condition, which must be satisfied by the initial 
data, is preserved by the evolution equations 33] . 



IV. SETUP FOR NUMERICAL 
COMPUTATIONS 

In this section we describe the setup we have used in 
the numerical simulations of the coupling of radial and 
axial non-radial oscillations of relativistic stars. A more 
detailed description of the numerical implementation can 
be found in |43j . 



A. Structure of the Numerical Codes 

It is clear from the characteristics of the perturbative 
framework described in Sec. [H] that the we have to im- 
plement the following computational steps, (i) Construc- 
tion, once and for all, of the TOV background; then, at 
each timestep: (ii) independent integration of the first- 
order radial and non-radial perturbative equations; (iii) 
updating, using the results of (ii), of the source terms in 
the perturbative equations describing the coupling; (iv) 
use the source computed in (iii) to integrate the equations 
for the coupling variables. In what follows we describe 
how these different steps are carried out. 

To solve the TOV equations (|11I12|) of hydrostatic 
equilibrium, we prescribe the value of the central density 
p c and integrate from the center to the surface, where 
the pressure vanishes. The integration is carried out by 
using a fourth-order Runge-Kutta method. 

Radial perturbations obey equations 119|) - (|23|) . It is 
well known (see, e.g. ^ij) that an accurate numerical 
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TABLE I: Convergence test for radial perturbations: cr*- 0,1 ' 
denotes the convergence rate in the L 2 norm. 



Variables: 


7 (o,i) 




5 (0,1) 


,,(0,1) 


CT (o,i) 


2.04 


2.06 


2.07 


1.53 



integration requires some attention at the stellar surface. 
Indeed, the low values that the speed of sound takes near 
the surface reduce the convergence rate of a second-order 
numerical scheme to first order. This problem can be 
solved by introducing, near the stellar surface, a tortoise 
fluid coordinate x |45| : 



dr = c c dx 



d r = (l/c s R 



(61) 



A uniform grid based on x (x-grid) will provide the neces- 
sary resolution near the surface by simulating a non uni- 
form mesh in the r coordinate. The transformation l|61(l 
has to be implemented in the TOV equations and in the 
equations for the radial perturbations (see Appendix 0. 

The are two ways to deal with the equations for the ra- 
dial perturbations. One is to use a free evolution scheme 
based on a purely hyperbolic formulation, which means 
to solve the equations l|C4|) - (|C8|l and to monitor errors 
through the Hamiltonian constraint l|C7|l . The second 
way is to use an elliptic-hyperbolic scheme, where we 
integrate Eqs. (|C4IC5|I for ijtbo) and 7( 1 '°) , and solve 
the Hamiltonian constraint l|C7l) for S^ 1 ' ) . In this sec- 
ond scheme the Hamiltonian constraint is satisfied at 
every time-step by construction. In both schemes, we 
need to integrate Eqs. (|C4IC5jl . for which we use a two- 
step Mac-Cormack algorithm with a predictor-corrector 
step to provide second-order accuracy both in space and 
time. When the metric perturbation S^ 1,0 ) is found from 
the evolution equation l|C8|l we use the same algorithm, 
whereas when it is found from the Hamiltonian con- 
straint i|C7|) we use a standard LU decomposition to solve 
the resulting tridiagonal linear system. Finally, the met- 
ric perturbation Tjvb ) is found, at every time step, by 
integrating equation (|C6|I by using the shooting method 
in order to implement the boundary condition 7y( 1,0 ) = 
at the surface. 

We use the second scheme for our computations, i.e. 
we integrate numerically equations (|C4I) - (|C7|) . The 
amplitude of the Hamiltonian constraint l|C7|l remains 
bounded by very small values for long-term evolutions 
and scales with the grid resolution as expected in a 
second-order convergent numerical evolution. The degree 
of convergence in the L 2 norm for the variables describing 
the radial perturbations are given in Table [I] 

Another test of the numerical code for the radial per- 
turbations consists in comparing the mode frequencies 
available in the literature, determined within a frequency 
domain approach, with those we obtain by applying a 
Fast Fourier Transformation (FFT) to the solution of our 
perturbative equations. In Fig. ^ we show the spectrum 
of the fluid velocity perturbation 'y( 1 '°) J where the radial 
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FIG. 1: Power spectrum of the radial fluid-velocity perturba- 
tion 7' 1 ' '. The upper panel displays simultaneously the FFT 
of eight different time evolutions. In each evolution a single 
radial mode has been excited. In the lower panel, radial oscil- 
lations are excited with an initial Gaussian pulse. The circles 
denote the frequencies of the radial modes determined from 
the eigenvalue problem (see Section TlV Bl , Note the different 
horizon scales used in the panels. 



oscillations have been excited with a selected eigenfrc- 
quency (upper panel) and with an initial Gaussian pulse 
(lower panel). The radial frequencies, excited during the 
time domain simulations, show an agreement to better 
than ~ 0.2 % with the eigenfrequencies determined in 
Section IIV Bl with a frequency domain code, and with 
the three known eigenfrequencies pEjj. 

In the case of the first-order axial perturbations equa- 
tions, as well as in the case of the non-linear coupling, 
near the stellar surface we don't have the problems ex- 
hibited by radial perturbations. Therefore, there is no 
need to introduce the tortoise coordinate and we can use 
an evenly spaced grid in r (r-grid) for the (0, 1) and (1, 1) 
perturbations. 

In order to build the source terms for the coupling be- 
tween first-order radial and non-radial perturbations, we 
need to interpolate the values of radial pulsations from 
the x-grid to the r-grid. We find a good agreement be- 
tween the evolved and interpolated radial variables when 
we use a x-grid whose resolution is twice that of the r- 
grid. 

We update in time the first-order axial perturbations 
by using equation H4fi[) for the axial master variable ^C' 1 ) 
and equation l)47|l for the axial velocity perturbation 
^(0,1) _ The integration of (|47|) is trivial: it yields a 
static profile for [3^°^ . For a nonzero velocity pertur- 
bation, the master equation (|46[) contains a static source 
term. Then, the solution of l|4(i[l can be written as: 



$(o,i) = $ 



(0,1) 



(0,1) 



, where 'I'j, ' 1 ^ 



is the solution of 



the homogeneous equation and describes the dynamical 
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degrees of freedom of the gravitational field, namely, it 
describes the gravitational wave content of the non-radial 
oscillations. 'J'p 0,1 '' is a particular solution of the full 
equations which can be chosen to be time independent, 
and is related to the dragging of the inertial frames due 



to the stellar rotation. The homogeneous part can be ob- 
tained numerically by using a standard Leapfrog scheme 
(see, e.g. H^). The particular solution 'J'p 0,1 ' satisfies the 
following ordinary differential equation: 



J 



p,rr 



„2A 



47rr(p — p) H ^- 



(0,1) 



_2A 



6M 1(1 + 1) 
4tt(p - p) + —3 — 2 — 



167rre 

r 



A 



(04) 



^(0,1) + e 2A ^4^+^^(0,1) 



= 0, 



(62) 



which can be discretized in space using a second-order 
Finite Differences approximation and written as a tridi- 
agonal linear system. Its solution can then be found by 
using a standard LU decomposition 0] . The particular 

solution ^p 0,1 " 1 is related to the gauge-invariant metric 

componen t fcp 0,1 ' 1 through the relation l|44|l . As it was 
argued in |37| . it describes the frame dragging of the in- 
ertial frames due to the presence of the axial velocity 
perturbation: 



where the equality fc 



(0,l)/m 



= h 



(o,i);m 



holds in the Regge- 



Wheeler gauge or when the field is time independent. 
The metric perturbation fcg ' 1 '''" 1 , as well as the frame 
dragging harmonic component to , can also be deter- 
mined directly from the following differential equation: 



-uj(r, 9)r 2 sin 2 I 



(0,1) 
St0 



(0,l)lm Q l 

J J. 



(63) 



,-2A 1,(0,1) 



47rr(p + p)fc 



(0,1) 



8tt(p- 



, AM 
P) + — - 



1(1 + 1) 



(64) 



r 



which follows from (|40|1 . (|44|l . It can be solved with the 
same method used for equation (|62J) . 

The linear non-radial perturbations are the combina- 
tion of a gravitational wave degree of freedom and of 
a static solution related to the dragging of the inertial 
frame. They correspond, respectively, to the dynamical 

Our 



solution ^^'^ and to the particular solution 'E r £ ' 1 -' 



tests show that and fc have the expected second 

order convergence, while 'J'p ' 1 ^ manifests a convergence 
of first order. This lower convergence rate is due to the 
discontinuity of ^^rr at the stellar surface. 

We have tested the performance of our numerical code 
for linear axial perturbations by comparing the results 
of our simulations of the scattering of an incident grav- 
itational wave onto a c omp act star with results in the 
literature (see, e.g. |47L l48lp. For the quadrupolar case 
(see Fig. 0), the scattered gravitational signal exhibits 
the characteristic excitation of the first iu-mode, followed 
by the ringing phase which is strongly damped by the 
emission of gravitational radiation. The behaviour of 



the signal at late times is known to be dominated by 
the gravitational-wave tails due to back-scattering of the 
waves by the spacetime curvature. The theoretically ex- 
pected time behaviour of the tails is vpC 3,1 ) ~ £-(2/+3) 
(see The behaviour we obtain from our simulations 

is shown in Fig. where by using a linear regression we 



obtain $ 



(0,1) 



t- 71 for I = 2 and *^ 0,1) 



t 



-9.26 



for 



1 = 3. We have checked that these exponents approach 
the theoretical values as we increase the evolution time. 
The spectral properties of the solution are extracted by 
Fourier transforming the signal. In Fig. |3 the power 
spectrum curves show the excitation of the lowest or- 
der w-modes. The frequencies associated with the peaks 
of these curves are in excellent agreement with the fre- 
quencies computed by the frequency-domain code used 
in 50] . Their values are: v w = 10.501 kHz for I = 2 and 
v w = 16.092 kHz for I = 3. The broad shape of the peaks 
is due to the short values of the w-modes damping time, 
t = 29.5 /is for I = 2 and r = 25.2 ps for I = 3. 

In order to test the numerical solutions of the equa- 
tion (|62|l for the axial master function \E , p°' 1 ' ) (for a par- 
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Log[t/(2M)] 

FIG. 2: Gravitational-wave signal determined by the scatter- 
ing of a Gaussian pulse on a TOV star. The exponentially 
damped sinusoidal oscillations corresponds to the I = 2 (top 
panel) and to the I = 3 (bottom panel) w-mode ringdown. The 
slope of the tail is in good agreement with the theoretical prc- 
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the only difference that the former contain source terms. 
The source terms E^ and E^ [their expressions are given 
in equations IjBlfl and l|B2|) respectively] are discretized 
by second-order centered finite difference approximations 
in the internal grid points and by second-order one-sided 
Finite Difference approximations at the origin and at the 
stellar surface. The computational domain for the ax- 
ial master equation l|58l) is the entire spacetime, but the 
sources Y,^ and E» have support only in the interior 
spacetime (inside the star). We have tested different nu- 
merical schemes to solve the equation for the master func- 
tion vl^ 1 ' 1 ). We have found that the scheme that provides 
the best accuracy depends on the physical setting that we 
are considering. In the case we are studying the coupling 
between the radial pulsations and differential rotation, 
our calculations are more accurate when we use an up- 
wind evolution algorithm for the interior and a leapfrog 
scheme for the exterior. In the case in which we study 
the scattering of an axial gravitational wave by a radi- 
ally oscillating star, the best results are obtained when 
we use a Leapfrog algorithm on the whole spacetime. The 
reasons for this difference are the different properties of 
the source term E^, and the junction conditions (|60fl in 
these two different physical scenarios (see below for a de- 
scription of the numerical implementation of the junction 
conditions). The evolution of the velocity perturbation 
pd' 1 ' is carried out by using an up-wind scheme (/3( 0,1 J 
does not need to be evolved since it is constant in time). 
The solution of the axial master equation H58[) and of the 
conservation equation l|59l) does not present any stabil- 
ity problem and the algorithm is first-order convergent. 
This is expected from the first order accurate numerical 
schemes used for the radially and differentially rotating 
configuration, and from the discontinuity of ^ at the 
stellar surface, due to the presence of the source terms 
only in the interior spacetime. 



FIG. 3: Spectrum of the first order axial perturbations due 
to the scattering of a Gaussian pulse on a TOV star. The 
peak correspond to the first ui-mode excitation, at frequencies 
v w = 10.5 kHz (I = 2) and v w = 16.092 kHz (1 = 3). 



ticular profile of the axial velocity perturbation described 
in subsection IIV B|) . we have compared the results with 
a different method. This alternative method consists in 
solving first the equation (|64[) for the metric perturba- 
tion k and then to obtain vE'p ' 1 ' 1 from the definition (|44|l . 
which for the stationary case becomes 

= ( 2 fc(°' 1) - rk^f) e-(*+ A > . (65) 

We have found that these two different computations 
agree to better than 2.3% (see Figure|SJ). 

The evolution equations for the non-linear coupling 
perturbations, equations (|58I59|I . have the same struc- 
ture as the equations for linear axial perturbations, with 



B. Construction of Initial Data 

In the following we describe how we prescribe initial 
data at the different perturbative orders: for the radial, 
the axial non-radial, and the coupling perturbations. 

The description of the radial perturbations requires to 
specify initial values for 7 ( 1 >°) , H (1 ^ , , and rf 1 ^ 

(although these variables are not independent since they 
are subject to the Hamiltonian constraint). In this paper 
we consider two different types of initial conditions: (i) 
We prescribe an initial profile for the fluid perturbations 
that corresponds to an eigenfunction associated with a 
single particular radial mode, (ii) We prescribe initially 
a Gaussian profile that will excite a broad range of normal 
radial modes. 

The way in which we implement both types of initial 
conditions consists of setting to zero the perturbations 
jj{i,o) ; ,5(1.0) ^ ^(1,0) an( j p rescr ibing a non-zero initial 
profile for jt 1 * ). For the first type of initial conditions 
(a single eigenmode) this implementation is equivalent to 
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the choice of a particular time origin. Indeed, as one can 
find out from equations (|19|) - (|21[1 and l|23|l . we can choose 
consistently an oscillation normal mode with eigenfre- 
qucncy a n to have the form 



l {1 ' 0) (t,r) 



7 (^)(r)cos(a„t), (66) 

i^MsinKi), (67) 
S£<°Hr)sin(* n t). (68) 



In this way we have that at t = only -yC 1 ' ) is non-zero. 
For the second type of initial data (Gaussian profile) our 
choice of implementation is not the most general one, it 
just represents a particular linear combination of eigen- 
modes. 

To construct the first type of initial data, we can set 
up the eigenvalue problem for the radial velocity pertur- 
bation 'y( 1,0 J , from the wave-type equation that can be 
derived from equations l(T§|l and (|2"U)l : 

(P(r)y^) r + Q(r)y^ - W(r)y% 0) = , (69) 



r [km] 




r [km] 



FIG. 4: Eigenfunctions for the radial velocity perturbation 
-y' 1 ' ' corresponding to the fundamental mode and the first 
seven overtones. 



where y&V = r 2 e -A 7 (i,o) = r^-*^ and ig the 
Lagrangian displacement of the radial oscillations. The 
functions W, P, Q are constructed from quantities associ- 
ated with the TOV background: 



r 2 W= {p+p)e 3A+ * . 
r 2 P={p + p)c 2 s e K+ ^ 



r 2 Q=(p + p) 



,r \ ,r 



8irpe 



2A 



(70) 
(71) 

3 A+3* ( ?2 ) 



By making the substitutions y(i'°) = y (r)e luJt and 
2(1,0) _ Py(y ^ in equation l|69|l we obtain the follow- 



ing pair of first-order equations 

,(1,0) = p-l z (l,0) _ 

-(u 2 W + Q)y^°K 



,(1,0) 



(73) 
(74) 



To complete the eigenvalue problem we need to prescribe 
boundary conditions at the origin and at the surface lo- 
cation. At the origin, from the regularity condition 
we must have the following behaviour 



(r) =y r 3 , (r) = zg' 0) 



(75) 



which combined with equation 1|73|) yields the following 
boundary condition: 



'0 ■ 



(76) 



The boundary condition at the surface comes from the 
vanishing of the Lagrangian pressure perturbation there. 
The resulting condition reads [see equation l|24|)]: 



(p+p)c 2 y^\ r=1 



0. 



(77) 



TABLE II: Eigenfrequencies (in kHz) of the first eight nor- 
mal modes of the radial perturbations. The second column 
corresponds to the results from calculation in the frequency 
domain (the solution of the Sturm-Liouville problem). The 
third column corresponds to the results from a time-domain 
calculation, after applying a FFT to the resulting time se- 
ries. The fourth column shows the first three normal modes 
computed by Kokkotas and Ruoff |2g|. 



Mode 


Frequency domain 


Time domain 


From [28] 


F 


2.138 


2.145 


2.141 


HI 


6.862 


6.867 


6.871 


H2 


10.302 


10.299 


10.319 


H3 


13.545 


13.590 




H4 


16.706 


16.737 




H5 


19.823 


19.813 




H6 


22.914 


22.889 




H7 


25.986 


25.964 





The eigenvalue problem can be solved numerically by 
using a standard relaxation method (for details see |46|L 
As discussed above, in order to obtain accurate evolu- 
tions of radial perturbations it is convenient to intro- 
duce the fluid tortoise coordinate (|61|l which provides 
the necessary resolution close to the stellar surface. The 
equations that one obtains for the eigenvalue problem in 
terms of the fluid tortoise coordinate have been given in 
Appendix |0 The eigenfrequencies and the associated 
eigenfunctions for the fundamental mode and the first 
seven overtones of the radial velocity 'y(i'O) are shown in 
table [H] and in Fig. 0] respectively. They have been nor- 
malized with respect to the absolute value of their maxi- 
mum. Our computations of the eigenfrequencies and the 
associated eigenfunctions have second-order convergence. 
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In the case of the first-order axial non-radial per- 
turbations, we need to prescribe initial data for vf^' 1 ) 
at its time derivative at an initial time t = t Q , say 

(\E , o°' 1 '\ <9 t \E , o°' 1 ' ) ) , and a certain profile for Z^ * 1 ) , or 
equivalently for [3^°^ (this quantity remains contant 
along the evolution). Specifying this profile is equiva- 
lent to prescribe a certain differential rotation law for the 
star [32ll37| . As we have discussed before, this law has an 
impact on the evolution of vf't 1 ' ) , and it induces a drag- 
ging of inertial frames. In this work we consider two types 
of initial data for the first-order axial perturbations: (i) 
we prescribe a certain rotation law Z^ ' 1 ) , and then we 
obtain a profile for $5 from the integration of equa- 
tion l|62|l which corresponds to the particular solution 
\I^, ' \ For the homogeneous part of vpf 1 ' 1 ), the dynami- 



cal one, we specify zero data: ^ , dt^ho') = (0,0) • 
(ii) We set the star rotation to zero: /3^ 0,1 ^ = 0. Then, we 
just set the particular solution ^p 0,1 ^ equal to zero and 
prescribe arbitrary profiles for the homogeneous part: 

In order to discuss the first type of initial data it is im- 
portant to mention that presently there is no sufficient 
knowledge on the laws of differential rotation for neutron 
stars. For nearly born neutron stars, one may obtain a 
description of the rotation law from the numerical sim- 
ulations of core collapse (see, e.g. [jHl E3- ^ n the 
Newtonian approach, a set of rotation laws has been in- 
troduced motivated by mathematical simplicity and the 
aim of satisfying Rayleigh's stability criterion for rotating 
inviscid fluids: 



d{g 2 Vt)/dg > 0. 



(78) 



of £1 only 54], that is, — In the slow rota- 

tion approximation this condition leads to the following 
expression: 



,■(0,1) 



t„.(o>i) 



-2$ 2 • 2 , 

e r sm I 



(fi-w) 



(80) 



which is a linear approximation in f2. The choice of the 
functional form of j(O), and hence of j'C ' 1 ) (Q), must sat- 
isfy the Rayleigh's stability criterion against axisymmet- 
ric disturbances for inviscid fluid: 



dl 

dn 



<o, 



where j is the specific angular momentum 



J=(p + P)^ 
Po 



(81) 



(82) 



and p is the rest mass density. The specific angular mo- 
mentum J is locally conserved during the axisymmetric 
collapse of a perfect fluid 0|- A common choice for j(SV) 
that satisfies these conditions 0, is the following: 



(83) 



where Q c is the angular velocity at the rotation axis and 
A is a constant that describes the amount of differential 
rotation. We can now find the rotation law from equa- 
tions {HDl and (113 • The result is: 

A z + e ZB r A sm 

In the Newtonian limit this equation reduces to the 7- 
constant rotation law used in Newtonian analysis 



where g = r sin 9 is the cylindrical radial coordinate and 
n = n(r,6) is the angular velocity measured by an ob- 
server at infinity and describes the stellar differential ro- 
tation. Among this family of Newtonian rotation laws 
there is one, the j-constant law, that has been incorpo- 
rated into a General Relativistic approach 0,11^1 by tak- 
ing into consideration the dragging of the inertial frames. 
In this work, we introduce a rotation law by prescrib- 
ing the velocity perturbation W^l from an expansion in 
(vector) spherical harmonics of the velocity perturbation 
of a slowly differentially rotating star. In the slow ro- 
tation approximation, the fluid velocity perturbation is 
given by 



fi(r,0) = 



A 2 n. 



A 2 



(85) 



A uniform rotating configuration with = f2 c is attained 
for high values of A (A — > 00). On the other hand, for 
small values of A, the law l|84|l describes, in the New- 
tonian limit, a configuration with constant angular mo- 
mentum. 

The only non-vanishing vector spherical harmonic 
components of ui ' 1 ^ are given by the expansion 



,(0,1) 



(86) 



,(0,1) 



o,< 



(0,1; 



(0,0,0,', 



(n-w)) , 



(79) 

where ft = Sl(r,d) is the angular velocity measured by 
an observer at infinity and describes the stellar differ- 
ential rotation, while the function u> = co(r, 9) describes 
the dragging of inertial frames associated with the stel- 
lar rotation. In the case of barotropic rotating stars, 
the integrability condition of the hydrostatic equilibrium 
equation requires that the specific angular momentum 
measured by the proper time of the matter is a function 



where Z?^' 1 ^ are the velocity perturbations (which de- 
pend on the coordinate r) that were introduced in equa- 
tion (|38|l [we have now restored the spherical harmonic 

indices (l,m)]. The different are obtained by using 

the properties of the inner product among the different 
elements of the basis of axial vector spherical harmonics. 
Their expression is 



0(0,1) = 
> J lm ~ 



1 



sm 



s 2 



(0,1) olm ab 
l a D b 7 • 



(87) 
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where 7 Q& is the contravariant metric on S . 

In order to determine the initial profile for the axial 
velocity perturbations we just have to introduce the ro- 
tation law H84[) into the form of the velocity perturba- 
tions We obtain 



,(o,i) _ 



t A 2 



A 2 



e -24> r 2 sm 2 , 



■ sin 2 e n r 



hn J 



where we have used equation I|tj3|) , which relates the met- 
ric perturbations with the frame dragging function lu . It 
is important to notice that the first term in (|88|) corre- 
sponds to the the Newtonian j-rotation law (up to the 
factor e — *), to which we will refer as the nearly New- 
tonian j-rotation law. The other terms account for the 
dragging of the inertial frames. To obtain the velocity 
perturbations 0^ from (|87|l we just need to introduce 
there expression Q88|) . In doing this, we can easily deter- 
mine the expression for the nearly Newtonian term, but 
we need to pay special attention to the relativistic cor- 
rections. The difficulty arises from the expansion in axial 
vector harmonics in l|88l) containing the unknowns, fc ™, 
of the differential equation (|64[) . Due to the form of (|88|l . 
the inner product in equation (|87|l will contain products 
of harmonics with different harmonic indices. In order 
to decouple these terms in the relativistic corrections we 
assume that the dominant contributions are provided by 
the metric perturbations k l Q m that have the same har- 
monic indices as in the nearly Newtonian rotation law. 
For the nearly Newtonian part, we have found solutions 
only for A > e~®( R ^ R s . From this relation and from the 



values assumed by the metric function <E> in the stellar 
model considered in this work (see subsection 1111 A|l we 
can give the following estimation for the allowed values 
of A: 



1 A 2 < 

2 < A 2 + e-^r 2 sin 2 9 ~ 



(89) 



Then, using these inequalities, the approximation that 
we obtain for the relativistic terms of equation l|88|l is 
given by: 



,(o,i) 



A 2 r 2 sin 2 fle- $ 
A 2 + e -2f r 2 sin 2 



a. 



1,1m qhn 



(90) 

where here a is a constant such that a £ (0.5, 1]. The 
coefficients of the expansion of this equation in axial vec- 
tor harmonics have the following form 



,(o,i) 

l <pM 



a e ^k'^S 1 ? for I even, 



0(0,1) 



for I odd , 



where the coefficients of the components with odd I are 
given by the following expression: 



0(0,1) 
Pio 



1) 



f l0 (x,A)+a e-*k l °, (91) 



where x = re (not to be confused with the fluid tor- 
toise coordinate). The functions f l0 (x,A) for I = 1 and 
I = 3 are given by the following expressions 
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(92) 
(93) 



which have been derived by imposing the condition A > 

R s e-*( R *l 

In the limit of A — ► 00, the functions / ;o behave as 
expected, that is, the nearly Newtonian part becomes: 



e-*ft c r 2 sin 2 6 for I = 1 , 



lim -4Jk=f lo sf 



A^OO y/l(l + l)' 







for I > 3 , 



which corresponds to a uniformly rotating configuration, 
where D, = Q c is the angular velocity as measured at in- 
finity. Since in this work we are mainly interested in the 
impact of the non-linear coupling of the stellar oscilla- 
tions in the gravitational wave emission, we will consider 



only the case I = 3, as the dipolar term I = 1 does not 
produce gravitational waves. 

In Fig. we plot the profiles for the axial velocity per- 



turbation p. 



(0,1) 

30 



It shows that the two solutions corre- 



sponding to the extreme values of the constant a dis- 
agree in an amount below 10%. In the following simula- 
tions we will use for simplicity only the nearly Newtonian 
term by obtaining results which are correct to better than 
10%. 

With regard to second choice of initial conditions for 
the axial non-radial perturbations, namely ft^ 0,1 ^ = , it 
is going to be used in this paper to investigate how the 
second-order metric perturbations corresponding to the 
coupling terms can be affected by the radial pulsations 
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FIG. 5: Profiles of the I = 3 component of the axial velocity 
perturbation /jJjo determined for a j-rotation law with A = 
10 km and a rotation period at the axis of T = 10 ms. The 
solid line corresponds to the axial velocity associated with the 
nearly Newtonian j'-law, while the dotted and dashed lines 
correspond to the velocity perturbations in 19111 with a = 0.5 
and a — 1 respectively. 



FIG. 6: The stationary axial master function 'J/p ' 1 ' relative 
to a nearly Newtonian j-rotation law with A = 15 km and 
a period T — 10 ms at the rotation axis (top panel). The 
solution of the equation 1621 is showed with the solid line 
while the dashed line denotes the solution found indirectly by 
first solving Eq. 1641 for the variable fco an d then using the 
definition 1651 . The frame dragging function ui 30 determined 
by the initial value 19111 is displayed in the bottom panel. 



of the star. Actually, it is well known that the axial 
components of the gravitational- wave signal emitted by a 
non-rotating compact star can only contain the imprint of 
the spacetime w-modes. The question we can ask is then 
how this particular signal looks like when it couples with 
the radial pulsations of the star. What we do in practice 
is to generate ui-modes at first-order and to study the 
corresponding signal at second-order (in the sector of the 
coupling of radial and non-radial modes) by looking at 
its dependence on the radial pulsations of the star. The 
w-mode oscillation is excited by the standard procedure 
of scattering a (narrow) Gaussian pulse of gravitational 
waves on the star. 

For the non-linear coupling perturbations we have to 
prescribe initial data for ijK 1 ' 1 ), * ( 4 M) and Z^ 1,1 -* at an 
initial time t — t Q . Since we are mostly interested in the 
non-linear effects generated from the coupling between 
linear perturbations and not in the evolution per se of 
second-order perturbations, we set these quantities ini- 
tially to zero. 

To sum up, combining the different types of initial data 
for the first-order perturbations we can study the follow- 
ing two different scenarios: (i) the scattering of a gravita- 
tional wave by a radially oscillating star, (ii) A radially 
pulsating star with with differential rotation. 

C. Boundary conditions 

The behaviour at the stellar center of the radial, ax- 
ial non-radial, and non-linear coupling perturbations is 
given in equations 1(25 (1 - 1(28 |) and 1)51(1. In order to deal 



numerically with the centrifugal term of the potential of 
the axial master equation, l(l-\-\)r~ 2 , we use a numerical 
grid where the first grid-point, r 1 , does not coincide with 
the center r = but instead is located at r 1 — Ar , being 
Ar the grid spacing. Then, we impose at r 1 the regu- 
larity conditions l(25 (l -H2H |) and l(51|) . For instance, in the 
case of the variable ^'( ' 1 ) j which behaves as ^ (f)r i+1 
[equation (|51|) ] . we assume that this behaviour is in a 
good approximation also valid at the grid points r 1 and 
r 2 , and then we determine the value of vpC' 1 ) a t r ^ f r0 m 
the following relation: 

*(°' 1 )(r 1 ) = *(°' 1 )(r 2 ) (^y +1 ■ (94) 

Let us consider now the boundary conditions at the stel- 
lar surface, where the pressure vanishes. In this way, 
the unperturbed stellar surface is E = {x | r(x) = R s } 
where R s is the stellar radius of the TOV model, where 
the background pressure p vanishes. Due to the ra- 
dial pulsations, the surface of the perturbed star does 
not coincide with the static surface, but is given by 
E = {x + A£ (1 ' 0) | x e E} , where is the La- 

grangian displacement vector of a fluid element. The 
axial non-radial modes only induce rotation of the star 
and hence, they do not change the shape of the stellar 
surface. Therefore, at first order we have to impose the 
vanishing of the total pressure on the perturbed surface 
E for the radial perturbations, which is equivalent to the 
relation l(24|) , and the continuity conditions on the static 
surface for the axial non-radial perturbations. 
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The radial perturbations are evolved using the x-grid 
(the one determined by the fluid tortoise coordinate), 
then we interpolate the solution to the r-grid (the one 
determined by the area radial coordinate). In terms of 
the fluid tortoise coordinate x, the surface condition 124|) 
is given in equation l(C14(l of AppendixEl with R x being 
the value of the fluid tortoise coordinate corresponding 
to the static stellar radius R s . This condition will be 
certainly satisfied by finite values for ^y^ 1 * ) and 7^ , 
since pressure, density and speed of sound go to zero 
at the surface for a polytropic equation of state. We 
then determine, at every time step, the finite value of 
jO-fi) on g U rface by using a second-order polynomial 
extrapolation. The surface boundary conditions for the 
other three radial variables are then directly determined 
from the perturbative equations. 

The solution of the axial master equation 1)46(1 is de- 
composed into two pieces: an homogeneous part ^j, ' 1 ^ 
that satisfies the homogeneous equation, and a static par- 
ticular solution i&p'^ that satisfies the equation (|62|1 . 
These equations are solved numerically on the whole spa- 
tial grid without imposing junction conditions. In this 
way, the continuity of and \E , p°' 1 ' ) and its first spa- 

tial derivatives, required by the junction conditions de- 
scribed in Section lill CI follows automatically. 

The junction conditions for the non-linear perturba- 
tions require some approximations (see |lOl |S| for a 
treatment of the junction conditions in a similar non- 
linear context). In an Eulerian gauge the perturbed sur- 
face will not coincide in general with the surface of the 
background equilibrium configuration. As a consequence, 
some perturbative quantities may take unphysical values 
near the surface during the expansions and contraction 
phases of the star. For instance, while the background 
enthalpy and density go to zero as we approach the static 
stellar surface, the associate radial perturbations have an 
oscillatory character with amplitudes that, near the sur- 
face, are bigger than the background values. Therefore, 
in a contraction phase, when the Eulerian radial pertur- 
bations of enthalpy and density are negative, the total 
mass-energy p + Sp^ 1 ' ^ may take negative values. Fur- 
thermore, the low densities which are present at the out- 
ermost layers of the star may also produce some numeri- 
cal errors in the simulations ■ These problems can be 
avoided if we do not impose the matching conditions for 
the non-linear perturbations at the static stellar surface, 
as we do for the linear perturbations, but at a hyper- 
surface that during the evolution remains always slightly 
inside the unperturbed star 0, |5tJ ■ In practice this pro- 
cedure is implemented as follows: (i) we estimate the 
amplitude of the surface movement due to radial pulsa- 
tions, (ii) We compute the linear (radial and non-radial) 
perturbations up to the static surface R s . (hi) We impose 
the junction conditions for the non-linear perturbations 
on a hypersurface that is always inside (even during the 
contraction phase) the unperturbed star, (iv) We es- 
timate the effects on the gravitational signal. In this 



sense, it is important to mention that in doing this ap- 
proximation we are removing the outer layers of the star, 
where the density is very low, which in practice means 
that we are neglecting less than one percent of the stel- 
lar mass, which does not produce significative changes 
in the waveforms and spectra of the gravitational signal. 
Moreover, by using an appropriate choice of initial con- 
ditions, the movement of the perturbed surface due to 
radial pulsations can be very small ~ 10~ 3 R S . There- 
fore, this approximative procedure leads to neglect only 
one or two grid points near the stellar surface. However, 
on the negative side, since the source term for the axial 
master equation l|58l) takes its maximum amplitude at 
the stellar surface, this removal of one or two grid points 
induces an error of about 5 % in the gravitational wave 
signal. 

The procedure we have just described allows us to sim- 
plify the junction condition IjtjUfl for the two initial con- 
figurations that we consider in this paper. In the case of 
the scattering of an axial gravitational wave by a radially 
pulsating star, the velocity perturbations /j^ ' 1 ) and p^ 1 ' 1 ' 
are set to zero. Therefore, by using the continuity of the 
first order perturbations, the junction condition 16(J() re- 
duces to the continuity of *b,r • We then evolve the axial 
master equation i|58[l for vpf 1 ' 1 ) on the whole numerical 
grid, with the source terms only present in the interior 
of the star. In this way the junction conditions on the 
matching surface are automatically satisfied. For the sec- 
ond situation that we study, the non-linear coupling in 
a radially pulsating and differentially rotating star, we 
can only use the continuity of the linear perturbations to 
reduce <|6(J[) to the following expression: 

e- A *^ 1 ' + l^rft 1 ^ . (95) 

In this case, the numerical procedure to solve the axial 
master equation i|58[) goes as follows: (i) In the interior 
spacetime we implement an up-wind scheme, (ii) At the 
matching hypersurface, the junction conditions provide 
the values of the axial function and its time and 

spatial derivatives, (iii) Finally, the values obtained at 
the matching hypersurface are used to evolve the axial 
master function tpC 1 * 1 ) in the exterior (the (1,1) Regge- 
Wheeler equation) by using a Leapfrog scheme. 



V. RESULTS 

The following section is devoted to discuss the results 
of time-domain numerical simulations describing the non- 
linear coupling of radial and axial non-radial oscillations 
of a static star in two different physical settings: (i) The 
scattering of a gravitational wave by a radially pulsating 
star ( subsection IV AJ) . (ii) A differentially rotating and 
radially oscillating star (subsection IV B|l . 
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FIG. 7: Waveforms (in logarithmic scale) from the scatter- 
ing of an I = 2 axial gravitational wave from a star oscillating 
radially in the fundamental mode (see Table HTl . The top 
panel presents the first-order axial non-radial master func- 
tion ^r' ' 1 ' (solid line) and the second-order coupling function 
ijft 1,1 ) (dashed line). The bottom panel exhibits the waveform 
* (0 ' i:i (solid line) together with the total signal * (0 ' 1) + 
(dashed line). 



A. Effects of radial pulsations on the scattering of 
a gravitational wave 

The aim of these simulations is to investigate the possi- 
ble signature of the radial oscillations of the star in scat- 
tered gravitational waves. The initial first order axial 
perturbation correspond to an I = 2 mode with vanish- 
ing velocity perturbation /?( 0,1 ). The profile of the axial 
master function ^(O' 1 ) is a Gaussian pulse 

* (0,1) (t o ,r) =A^e-^ r - r ^\ (96) 

centered at r$ = 20 km with amplitude A^ ' 1 ) = 0.1 and 
width parameter q — 1.25. The star oscillates in one of 
the radial eigenmodes corresponding to the frequencies 
of Table HU 

In Figs.[7|and|Sl we present the I — 2 waveforms associ- 
ated with the first-order axial non-radial master function 
^r(o.i) anc j t ne mas ter function describing the cou- 

pling. We find that the correction to the signal coming 
from the coupling is less than 2 % when the radial pulsa- 
tions are excited by the F-modes, and less than 0.1 % for 
higher overtones. These corrections do not modify the 
properties of the waveforms and spectra associated with 
the first-order perturbations. We find similar results for 
higher overtones and their combinations. In summary, 
in this particular scenario we have not found any signifi- 
cant amplification of the gravitational wave signal due to 
the coupling of the radial and axial non-radial first-order 
perturbations. 




Log|t/(2M)| 



FIG. 8: Waveforms (in logarithmic scale) from the scatter- 
ing of an / = 2 axial gravitational wave from a star oscil- 
lating radially in the first overtone HI (see Table IHl . The 
top panel presents he first-order axial non-radial master func- 
tion i]/' ' 1 ) (solid line) and the second-order coupling function 
vjH 1 ' 1 ) (dashed line) . The bottom panel exhibits the waveform 
't> {0 ' 1) (solid line) together with the total signal * (0 ' 1) + 
(dashed line). 



B. Coupling between radial pulsations and 
differential rotation 

The aim of the simulations we present in this subsec- 
tion is to investigate the effect of the coupling of radial os- 
cillations with differential rotation, as described by axial 
non-radial first-order perturbations, in the gravitational 
wave signal. The axial rotation of the fluid is described 
with a good approximation by the nearly Newtonian j- 
constant rotation law described in the subsection IIV Bl 
All the simulations in this subsection have a first-order 
rate of convergence and are long term stable. 



1. Single radial mode oscillation 

We start with simulations in which we consider single 
mode excitation of the radial oscillations, in particular 
the first overtone HI. The initial velocity perturbation 
ry(i)0) i s taken to be [see equation JSS}] 

7 (i,o) = G (i,o) 7 (i,o) t (97) 

with G^' 0) = 0.001. For this amplitude we can be sure 
that the region of the spacetime that the motion of the 
surface covers is confined to a narrow slice around the 
static equilibrium configuration. Therefore, the issues 
related to negative values of the mass-energy density in 
an Eulerian description can be solved approximatively as 
we explained above. 
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FIG. 9: Coupling between radial mode oscillations (overtone 
HI) and differential rotation (with rotation period T = 10 
ms). Waveform vl/' 1 ' 1 ' for / = 3 multipole (top panel). A 
(transient) excitation of the first I — 3 ui-mode at early times 
is present, followed then by a periodic oscillation driven by the 
radial mode. The Fourier spectrum of the periodic oscillation 
(bottom panel) is dominated by the frequency i/hi — 6.8 kHz 
corresponding to the HI radial pulsation. 



The values of the differential rotation parameter A 
and the angular velocity at the rotation axis in these 
simulations are A — 15 km and Vl c = 2.09 x 1CP 3 
km . This latter (corresponding tp a period of 10 ms) is 
small compared with the mass shedding limit value Sl K : 
Q, C /Q, K = 6.45 x 10~ 2 . The profile of the axial velocity 
perturbation /J^ 0,1 ) for I = 3 is given in equation (|93|) and 
depicted in Fig. [5] For this initial configuration of the dif- 
ferential rotation, the coupling non-linear perturbations 
are dominated by the harmonic I — 3 . 

In Fig. we show the I — 3 waveform (top panel) ex- 
tracted at 100 km apart from the stellar center: after 
an initial transient phase dominated by w-modc ring- 
ing (see below), the vlK 1 ' 1 ) waveform exhibits a monocro- 
matic oscillation at the frequency of the HI mode; i.e., 
~ 6.8 kHz, as confirmed by the Fourier spectrum of 
the signal (bottom panel). The presence of the ui-mode 
in the initial phase of the waveform is related to the 
choice of initial data that we have done for simplicity; 
i.e., = # ( t M) = = 0. This is, in general, non 

consistent, because these three quantities are not inde- 
pendent; indeed, a correct set up of the initial data for 
v]/(M) would require the solution of the "coupling" axial 
constraints (not reported here) in a way similar to what 
is usually done for the first-order polar and axial pertur- 
bations (see for example Ref. |42|). If this is not done, the 
system requires to relax to the "correct" solution through 
the release of a GWs burst of w-modes. However, since 
our interest is focused here on the stationary phase driven 
by the radial oscillations of the star, we decided to avoid 
the complications related to correctly capture the initial 



transient. 



2. Amplification effects 

In this section we extend the previous analysis to simu- 
lations where the radial pulsations are excited using com- 
binations of the fundamental mode and higher overtones. 
In general, the waveforms and spectra present the same 
features as the HI case described above. However, an in- 
teresting amplification is observed in simulations where 
the radial oscillations contain a frequency close to the 
one of the axial w-mode. In Fig. ^| we show the axial 
master function ip' 1 ' 1 ) for six runs using different radial 
overtones. It is clear there that vpC 1 . 1 ) increases in ampli- 
tude as we increase the order of the radial normal mode 
from the first to the fourth overtone, while it decreases 
as we go to the fifth and sixth overtones. The ampli- 
tude of vfrC 1 ' 1 ) for the H4 case is about sixteen times that 
of vfK 1 ) 1 ) for the HI case. For this stellar model, the 
/ = 3 spacctime mode has frequency 16.092 kHz which 
is between the frequencies of the third and fourth ra- 
dial overtones, 13.545 kHz and 16.706 kHz respectively. 
Moreover, this effect takes place despite the energy and 
the maximum displacement of the surface of the radial 
modes decrease proportionally to the order of the radial 
modes (see Table IIII|1 . The interpretation of this am- 
plification as a sort of resonance effect is supported by 
the structure of the axial master equation (|58|l and its 
analogy with a forced oscillator. The amplification arises 
when one of the natural frequencies, determined by the 
form of the axial potential V [see equations (|50l I58|) ]. 
is sufficiently close to the frequency (or frequencies) as- 
sociated with the external force. In the case when the 
non-radial perturbations are static and just describe dif- 
ferential rotation, the frequencies corresponding to the 
external force are determined by the structure of the ra- 
dial oscillations. This picture is confirmed by the fact 
that this amplification effect disappears when the axial 
potential is arbitrarily removed. In addition, the fluid 
velocity perturbation (3^ x,1 \ which satisfies a completely 
different equation (|59|l , does not show any such amplifica- 
tion and instead decreases proportionally with the order 
of the radial mode (see Fig.lTTI). 

In agreement with these results, the power emitted in 
gravitational waves to infinity exhibits the same amplifi- 
cation too (see Fig. 112(1 . There are obviously two differ- 
ent contributions to the power, one due to the w-mode 
excitation related to the initial data choice and the sec- 
ond related to the periodic emission driven by the stellar 
radial pulsations. The relative strength of these two con- 
tributions to the power changes as we change the radial 
eigenmode. The periodic emission increases for the first 
eigenmodes and reaches its maximum for the H4 over- 
tone; it decreases then for H5 and H6. 

Let us now look at the coupling between differential 
rotation and radial pulsations in the fundamental mode 
F. The vJH 1 '!) waveform is reported in Fig. ED After the 
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TABLE III: Quantities associated with radial normal modes and their coupling to the first-order axial differential rotation: 
Energy, E^'°\ and maximum stellar surface displacement of the radial eigenmodes [From initial conditions 19711 with 

G (1 ' 0) = 0.001]; power, .Egjj , emitted in gravitational waves to infinity from the coupling between the radial eigenmode and 
the axial differential rotation; estimated values of the damping times, r!^' 1 ^ ; and number of oscillation periods, N osc , that takes 
for the non-linear oscillations to radiate the total energy initially contained in the radial modes. 



Normal Mode: 


F 


HI 


H2 


H3 


H4 


H5 


H6 


E£' 0) [10- 8 km] 

d 1,0) H 

<Eg^> [ID" 13 ] 

T 3 y [ms] 

N 

J * OSC 


35.9 

12.65 
1.85 x 10 -6 
6.49 x 10 9 
1.39 x 10 10 


4.2 
4.02 

6.83 x 10 -2 
20.49 x 10 3 
1.408 x 10 5 


1.37 
2.66 
4.85 
94.15 
971.58 


0.62 
2.02 
56.03 
3.69 
49.99 


0.34 
1.64 
157.01 
0.72 
12.07 


0.21 
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19.08 
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0.14 
1.19 
4.45 
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240.77 
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FIG. 10: Comparison of six vj/' 1,1 ' waveforms for the I — 3 
multipole. The differential rotation law is the same as the one 
used in Fig. [5] The radial pulsations considered correspond to 
single mode oscillation from HI to H6 overtone. These plots 
show that a resonance effect take place in vp^ 1 ' 1 ). See text for 
a discussion. 
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FIG. 11: Time evolutions (the same runs as in Fig. of 
the second-order fluid velocity perturbation /3*- 1 ' 1 ' . The per- 
turbation has been averaged over the star at each time 
step and plotted here. The dashed line in the left top panel 
(HI) reports the perturbation for the H2 eigenmode and is 
superposed for comparison. 



initial w-mode transient, the signal is dominated by a 
periodic pulsation at the radial F-mode frequency. We 
notice that the amplitude of the oscillation is smaller 
than in the previous cases due probably to the fact that 
the radial F-mode frequency is not enough close to the 
first w-mode frequency. The periodic part of the signal 
exhibits high frequency oscillations associated with the 
low density region near the stellar surface: in Fig. 1131 
we depict (with a dashed line) the waveform obtained 
when the junction conditions have been imposed at a 
radial coordinate r — 8.64 km, (which is equivalent to 
neglect 0.3 % of the stellar mass) and with a solid line the 
case when the matching surface stands at r = 7.75 km 
(corresponding to neglecting a 6.3% of the total mass). 
In the later case the numerical noise is reduced and the 



waveform looks closer to the expected form. 

Our perturbative approach does not include backreac- 
tion (see, e.g. [HMIE3)i m other words, it does not account 
for the damping of the radial oscillations or the slowing 
down of the stellar rotation due to contribution of the 
non-linear coupling to the energy and angular momen- 
tum loss in gravitational waves. Backreaction could be 
studied by looking at higher perturbative orders. Never- 
theless, we can provide a rough estimate of the damping 
time of the radial pulsations by assuming that the energy 
emitted is completely supplied by the first-order radial 
oscillations, and that the power radiated in gravitational 
waves is constant in time. In this way, the damping time 
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FIG. 12: Emitted power in gravitational waves pv 1 , 1 ) 



p(M) 



for the six runs of Fig. 1101 



is given by the following expression: 



r (i,i) 

'im 



E, 



(1,0) 



<^> 



where -En 1,0 ' is the energy of a radial eigenmode (see Ta- 
ble llllf) . and < > is the averaged value of the non- 

linear coupling contribution to the power emitted. The 
results for t^' 1 ^ are shown in Table ITTT1 Moreover, in the 
last row ofTablelnllwe give an estimation of the damping 
of the radial pulsations associated with a certain radial 
eigenmode in terms of the number of oscillation cycles: 



r (l,i) 

'lm 

T 



(99) 



where T n — v^ 1 , with v n being the eigenfrequency of the 
radial eigenmode. 

Finally, it is important to make some comments on the 
difference between the two scenarios we have studied: (i) 
scattering of gravitational waves by a radially oscillating 
star, where we did not find any amplification or reso- 
nant effect, and (ii) the coupling of differential rotation 
and radial pulsations, where we have observed amplifi- 
cation when the frequencies associated with these two 
different first order perturbations are close. This differ- 
ence is essentially due to the different properties of the 
source terms in these two cases, which are determined 
by the first order radial and axial non-radial perturba- 
tions. In the first scenario, the source acts typically for 
a relatively short period of time, which is given by the 
travel time of a gravitational wave across the star. Af- 
ter the gravitational waves have been scattered, the first 
order signal still present inside the star decreases accord- 
ing to the power law for gravitational wave tails. As a 
consequence, there is a ringing phase in the second order 




FIG. 13: Waveforms generated by the coupling between axial 
differential rotation and the radial fundamental mode. We 
show superposed the waveform with the matching condition 
at r = 8.64 km (dashed line) and at r = 7.75 km (solid line). 
See text for a discussion. 



perturbations but it is well below the amplitude of first 
order signal. On the other hand, in the second scenario, 
the source terms act periodically on the star forever, as 
this model does not contain the back-reaction. As a re- 
sult, the source has more time to couple with the axial 
non-radial perturbations. 



VI. CONCLUSIONS AND DISCUSSION 

In this paper we have investigated non-linear effects 
in the dynamics of compact stars due to the coupling of 
radial and axial non-radial oscillations, and the potential 
impact that this mechanism may have in the emission 
of gravitational radiation. To that end, we have used 
and extended the multi-parameter non-linear perturba- 
tive scheme introduced in 0, [2^ , and further developed 
in [23( for the case in which the coupling is between radial 
and polar non-radial oscillations. 

The non-linear perturbative equations and the gauge- 
invariant character of the perturbative quantities de- 
scribing the coupling have been derived by using the 2- 
parameter perturbation theory in connection with the 
GSGM formalism [Hill HI- As expected, in the stellar 
interior, the perturbative variables describing the cou- 
pling satisfy inhomogeneous linear equations. The form 
of these equations is such that the associated homoge- 
neous equations are constructed from the same linear op- 
erators as the equations for the first order axial non-radial 
perturbations. The source terms are quadratic in the 
first-order radial and axial non-radial perturbative quan- 
tities. In the exterior we do not have source terms and 
hence, the dynamics is described by Regge- Wheeler-type 
equations. Interior and exterior communicate through 
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the junction conditions at the surface of the star. We 
have discussed in detail the initial data setup used in the 
simulations presented in this paper as well as the different 
boundary conditions. The initial data for the radial os- 
cillations corresponds to specific eigenmodes. In the case 
of the axial non-radial perturbations we have used two 
types of initial data: (i) static profiles representing differ- 
ential rotation of the star from the relativistic j-constant 
rotation law, and (ii) Gaussian packets containing a num- 
ber of axial eigenmodes. Regarding the boundary condi- 
tions, of particular importance are the conditions at the 
stellar surface. In order to avoid negative values of the 
mass energy density near the surface due to the Eulerian 
description of the radial and non-radial polar perturba- 
tions, we have adopted an approximation already used 
in the literature |10L |4J, |57| , that consists in removing 
outer layers of the neutron star, which is equivalent, in 
most of the cases we deal with, to neglect less than one 
percent of the total gravitational mass of the star. This 
approximation leads to a description of the second-order 
coupling effects which is accurate to better than five per 
cent. 

We have presented a computational framework for in- 
vestigating, in the time domain, the evolution of the axial 
non-linear perturbations describing the coupling. Our 
numerical codes are based on finite difference methods 
and standard explicit evolution algorithms. Their struc- 
ture is hierarchical: First, we solve the TOV equations; 
then, we evolve a time step the first-order radial and non- 
radial perturbations, and with the result we update the 
source terms for the non-linear perturbations; eventually, 
we evolve them in the same time step. 

We have then studied two different physical scenarios: 
(i) The scattering of a gravitational wave by a radially 
pulsating star, (ii) A differentially rotating and radially 
oscillating star. The simulations for the first scenario 
have shown that the correction onto the gravitational 
wave signal due to the non-linear coupling are of the 
order of 2 % or less when the radial oscillations are in 
the fundamental mode, and of the order of 0.1 % or less 
when we consider higher radial overtones. Moreover, the 
properties of the waveforms and spectra are essentially 
unchanged with respect to the first-order ones. The sim- 
ulations for the second scenario have provided more in- 
teresting results. The waveforms we obtain have the fol- 
lowing pattern: at the early stages of the evolution they 
present an initial w-mode burst related to the initial data 
choice. At later stages, the signal becomes periodic and 
is driven by the radial pulsations through the sources. 
This is confirmed by an analysis of the spectra, which 
show that the gravitational wave signal contains the fre- 
quencies of the radial eigenmodes prescribed in the ini- 
tial data. More interestingly, we have found a resonance 
effect that takes place as the frequencies of the radial 
eigenmodes considered get close to the frequency of the 
first lu-mode. For the stellar model used in this paper, 
the amplitude of the gravitational wave signal associated 
with the fourth radial overtone (the one with the closest 



frequency to the first u>-mode) is about three orders of 
magnitude higher than that associated with the radial 
fundamental mode. We have also given a rough estimate 
of the damping times of the radial pulsations due to the 
gravitational wave emission. The values found depend 
strongly on the presence of resonances. In the case of 
differential rotation with a 10 ms rotation period at the 
axis and a differential rotation parameter A — 15 km, 
the fundamental radial mode gets damped after about 
ten billion oscillation periods, while the fourth overtone 
after only ten. This is not surprising, and shows that 
the coupling near resonances is a very effective mecha- 
nism for extracting energy from the radial oscillations. 
In this sense, it is important to mention that the detec- 
tion of such a gravitational wave signature would provide 
important information about the stellar physical proper- 
ties, specially because it has information on radial pulsa- 
tions, which can be determined easily for a large class of 
equations of state. 

The work presented in this paper represents a first step 
in the study of the non-linear coupling of oscillations of 
compact relativistic stars and its impact for gravitational 
wave physics. The next step in this study is the descrip- 
tion of the coupling between radial and polar non-radial 
oscillations. The spectrum of the polar non-radial pertur- 
bations is richer than the axial case due to the presence of 
the fluid modes, which may have frequencies lower than 
the spacetime modes and a longer gravitational damping. 
As a result, we may expect a more effective coupling with 
the radial pulsation modes and more channels for possible 
resonant amplifications. Other extensions of this work 
include the consideration of more realistic models of the 
stellar structure, by taking into account the effects of ro- 
tation, composition gradients, magnetic fields, dissipative 
effects, etc. In particular, it would be interesting to inves- 
tigate a proto-neutron star model to compare the damp- 
ing rates due to gravitational radiation produced by the 
coupling of oscillations with the strong damping induced 
by the presence of a high-entropy envelope surrounding 
the newly created neutron star. By including rotation 
we may find new interesting non-linear effects due to the 
different behaviour of radial and non-radial modes in a 
rotating configuration. While the radial modes are only 
weakly affected so that their spectrum is essentially the 
same as in a non-rotating configuration (when scaled by 
the central density), the non-axisymmetric modes man- 
ifest a splitting similar to the Zeeman effect in atomic 
spectra. The rotation removes the mode degeneracy in 
the azimuthal harmonic number of the non-rotating case. 
The details of this frequency separation depend on the 
stellar compactness and rotation rate. It may then hap- 
pen that for a given stellar rotation rate and compactness 
the non-radial frequencies cross the sequence of radial fre- 
quencies ^3 so that the frequencies of these two kinds 
of modes are comparable, and possible resonances or in- 
stabilities could influence the spectrum and the gravita- 
tional wave signal. Finally, as we have already mentioned 
above, it would be also interesting to explore the impact 
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of gravitational back-reaction on the dynamical coupling 
mechanism (see, e.g. [S|||5|j). 



where V^ 1,1 ^ are scalar functions on M 2 . Then, the met- 
ric perturbations h^' 1 ^ and h^ 1 ' 1 ^ , under the gauge trans- 
formation IjAlfl . transform as follows: 



;(-r,i) _ 
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APPENDIX A: GAUGE INVARIANCE 

In this section we show the gauge-invariant character of 
the (1,1) perturbative quantities H55fl - I|57[l (see for more 
details |23j). Gauge transformations and gauge invari- 
ance in 2-parameter perturbation theory have been stud- 
ied in pll I22] ] . The gauge transformation of a perturba- 
tion of order (0, 1) of a generic tensorial quantity T is 
given by 



f (0,1) = y(0,l) 



£ T (o,o) 

S (0,D 



(Al) 



whereas a perturbation of order (1,1) of T transforms 
according to (see p2|l 



j-(i.i) 



+ £ e 



+ Ut +{£* }) r im , (A2) 



where { , } stands for the anti-commutator {a, b} = a b + 
ba. In our work we fix the gauge for the radial (1,0) 
perturbations. Then, the transformation l|A2(l becomes 



2 T/ (/,l) 



(A5) 



~ h (I,l) = h (I,l) +r 2y(I,l)_ (A6) 

and the (1, 1) energy-momentum tensor perturbations as 



5? 



1,1) 



- 2 \pVW + pclu^vW 



(A7) 
(A8) 

it is 



Then, by using the transformation rules (|A5I 
straightforward to prove that the variables k^' 1 ^ , L^' 1 ^ , 
and ZA 1 ' 1 ) are gauge-invariant in the sense explained 
above. One can also apply the same procedure to study 
the gauge invariance of the coupling perturbation vIa 1 ' 1 ) . 
By using the ansatz l|54|) in combination with its expres- 
sion i|4T)|) at the (0, 1) order we get 



-k, 



^(1,0)^(0,1) 



(1,1) 



2k, 



(1,1) 



,-(*+A) 



(A9) 



Since vjK 1 ' 1 ) depends linearly on the gauge-invariant 
quantity and contains the product of the radial 

perturbation yy^ 1 ' ) with the first-order gauge-invariant 
axial perturbation v^ , 1 ) , we have that, for a fixed ra- 
dial gauge, the variable 'jK 1 ' 1 ) is gauge-invariant, in the 
sense that it is invariant under the transformations (|A1|) 
and 1A3II. 



Finally, the gauge invariant character of the velocity 
perturbation /J^ 1 ' 1 ) follows from the fact that u a — and 

that Ua = . 

? (0,1) 



<f(M) = 7-(i 



? (0,1) s (i,i) 



(A3) 



To start with, let us expand in odd-parity vector har- 
monics the generators ^ and ^ ^ of the gauge trans- 
formations: 



I = 0,1. 



(A4) 



APPENDIX B: SOURCE TERMS FOR THE (1,1) 
PERTURBATION EQUATIONS 



In this Section we give the expressions of source terms 
and So that appear in the equations 1)58 (1 - 1|59 |) for the 
(1, 1) perturbations: 
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APPENDIX C: INTRODUCING THE FLUID 
TORTOISE COORDINATE 

Here we provide the expressions for the TOV equations 
and the equations for the radial perturbations when we 
use the fluid tortoise coordinate x [see Eq. (|fc> II) ]. The 
TOV equations (|f flf2f) become: 



P,x = ~(P + P)®,x 

M- = 4nr 2 pc s , 



(C2) 
(C3) 



The evolution equations (jf 9|l - 12f|l for the perturbations 
ij(i>°J , 7 (!'°) and S^ 1,0 ) , the equation for the metric per- 
turbation r^ 1,0 ) (122[| . and the Hamiltonian constraint 123() 
take the following form: 
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The equations of the eigenvalue problem for radial per- 
turbations read as follows: 



= c s P~^ , (09) 
z%<V = ~c s (u 2 W + Q)y^\ (CIO) 



where i/ 1 ' ) = r 2 e A ^(i>°) = 
functions W, P, Q are given by: 

r 2 W= (p + p)e 3A+ *, 
r 2 P={p + p)c 2 e A+ ^ , 

r 2 Q = (p + p) 1 ■ 



r 2 e *£ r t , and where the 
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The boundary conditions transform at the origin and on 
the surface as follows: 



Vo 



{p + P)c s e 



S 3P 
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0. (C14) 
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